#ifndef _LBWRAPPER_HH_
#define _LBWRAPPER_HH_

#include "olb3D.h"
#ifndef OLB_PRECOMPILED // Unless precompiled version is used,
  #include "olb3D.hh"   // include full template code
#endif
#include "../sMathLib/VectorX/VectorX.hpp"
#include "BoundaryType.h"
#include "LBMUnit.hpp"
#include "StringUtility.h"

#include <string>
#include <sstream>
#include <vector>

using namespace olb;
using namespace olb::descriptors;
using namespace olb::graphics;

///////////////////////// Implement by liao /////////////////////////
template<typename T, template<typename U> class Lattice>
struct BackupCellF : public WriteCellFunctional<T,Lattice>
{
    virtual void apply(Cell<T,Lattice>& cell) const
	{
		cell.backupF(); 
	}
};

template<typename T, template<typename U> class Lattice>
struct TurnToBackupCellF : public WriteCellFunctional<T,Lattice>
{
    virtual void apply(Cell<T,Lattice>& cell) const
	{
		cell.turnToBackupF(); 
	}
};

template<typename T, template<typename U> class Lattice>
struct InitExternalForce : public WriteCellFunctional<T,Lattice>
{
	InitExternalForce(T extF[3])
	{
		for(int dim=0; dim<3; dim++)
			_extF[dim] = extF[dim]; 
	}
    virtual void apply(Cell<T,Lattice>& cell) const
	{	
		T* force = cell.getExternal(Lattice<T>::ExternalField::forceBeginsAt); 
		for(int dim=0; dim<3; dim++)	
			force[dim] = _extF[dim];
	}
	T _extF[3];
};
///////////////////////// Implement by liao /////////////////////////

template<typename T, template<typename U> class Lattice>
void continueGeometry
(
	BlockStructure3D<T, Lattice>& lattice,
	sMathLib::VectorX<int, 3> domain,
	T omega,
	Dynamics<T, Lattice>& bulkDynamics,
	OnLatticeBoundaryCondition3D<T, Lattice>& bc,
	BoundaryInput<T>& bInput
)
{
	std::cout << "LBWrapper.hh ver. 2011.5.18.1" << std::endl;

	const int nx = lattice.getNx();
	const int ny = lattice.getNy();
	const int nz = lattice.getNz();

	std::cout << nx << " " << ny << " " << nz << std::endl;

	bool modifyBC = false;

	lattice.defineDynamics(0,nx-1, 0,ny-1, 0,nz-1, &bulkDynamics);

 	if(bInput.getBoundary("Z-Left-hand-side").first  == "Pressure-Boundary")
 	{
 		bc.addPressureBoundary2N(   1,nx-2,   1,ny-2,   0,   0, omega); //set pressure B.C. on bottom
 		std::cout << "Accept the Pressure Boundary at 2N. " << std::endl;
 		modifyBC = true;
 	}
 	else if(bInput.getBoundary("Z-Left-hand-side").first  == "Velocity-Boundary")
 	{
 		bc.addVelocityBoundary2N(   1,nx-2,   1,ny-2,   0,   0, omega);
 		std::cout << "Accept the Velocity Boundary at 2N. " << std::endl;
 		modifyBC = true;
 	}

 	if(bInput.getBoundary("Z-Right-hand-side").first  == "Pressure-Boundary")
 	{
 		bc.addPressureBoundary2P(   1,nx-2,   1,ny-2,nz-1,nz-1, omega); //set pressure B.C. on top
 		std::cout << "Accept the Pressure Boundary at 2P. " << std::endl;
 		modifyBC = true;
 	}
 	else if(bInput.getBoundary("Z-Right-hand-side").first  == "Velocity-Boundary")
 	{
 		bc.addVelocityBoundary2P(   1,nx-2,   1,ny-2,nz-1,nz-1, omega);
 		std::cout << "Accept the Velocity Boundary at 2P. " << std::endl;
 		modifyBC = true;
 	}

 	//bc.addVelocityBoundary2N(   1,nx-2,   1,ny-2,   0,   0, omega);
 	//bc.addVelocityBoundary2P(   1,nx-2,   1,ny-2,nz-1,nz-1, omega);

 	bc.addVelocityBoundary0N(   0,   0,   1,ny-2,   1,nz-2, omega);
    bc.addVelocityBoundary0P(nx-1,nx-1,   1,ny-2,   1,nz-2, omega);
    bc.addVelocityBoundary1N(   1,nx-2,   0,   0,   1,nz-2, omega);
    bc.addVelocityBoundary1P(   1,nx-2,ny-1,ny-1,   1,nz-2, omega);

    bc.addExternalVelocityEdge0NN(   1,nx-2,   0,   0,   0,   0, omega);
    bc.addExternalVelocityEdge0NP(   1,nx-2,   0,   0,nz-1,nz-1, omega);
    bc.addExternalVelocityEdge0PN(   1,nx-2,ny-1,ny-1,   0,   0, omega);
    bc.addExternalVelocityEdge0PP(   1,nx-2,ny-1,ny-1,nz-1,nz-1, omega);

    bc.addExternalVelocityEdge1NN(   0,   0,   1,ny-2,   0,   0, omega);
    bc.addExternalVelocityEdge1NP(nx-1,nx-1,   1,ny-2,   0,   0, omega);
    bc.addExternalVelocityEdge1PN(   0,   0,   1,ny-2,nz-1,nz-1, omega);
    bc.addExternalVelocityEdge1PP(nx-1,nx-1,   1,ny-2,nz-1,nz-1, omega);

    bc.addExternalVelocityEdge2NN(   0,   0,   0,   0,   1,nz-2, omega);
    bc.addExternalVelocityEdge2NP(   0,   0,ny-1,ny-1,   1,nz-2, omega);
    bc.addExternalVelocityEdge2PN(nx-1,nx-1,   0,   0,   1,nz-2, omega);
    bc.addExternalVelocityEdge2PP(nx-1,nx-1,ny-1,ny-1,   1,nz-2, omega);

    bc.addExternalVelocityCornerNNN(   0,   0,   0, omega);
    bc.addExternalVelocityCornerNNP(   0,   0,nz-1, omega);
    bc.addExternalVelocityCornerNPN(   0,ny-1,   0, omega);
    bc.addExternalVelocityCornerNPP(   0,ny-1,nz-1, omega);
    bc.addExternalVelocityCornerPNN(nx-1,   0,   0, omega);
    bc.addExternalVelocityCornerPNP(nx-1,   0,nz-1, omega);
    bc.addExternalVelocityCornerPPN(nx-1,ny-1,   0, omega);
    bc.addExternalVelocityCornerPPP(nx-1,ny-1,nz-1, omega);

 	for (int iX=0; iX != nx; ++iX)
 	{
 		for (int iY=0; iY != ny; ++iY)
		{
 			for (int iZ=0; iZ != nz; ++iZ)
			{
				T vel[] = { 0., 0., 0. };
				lattice.get(iX,iY,iZ).defineRhoU(1., vel);
				lattice.get(iX,iY,iZ).iniEquilibrium(1., vel);
			}
		}
 	}
/////////////////////// Initialize Field ///////////////////////

/////////////////////// Restore Field ///////////////////////

 	std::cout << " CP = " << 0.1 << std::endl;

 	using namespace std;

 	vector< vector< vector < vector< T > > > > _field3d;

 	std::cout << " CP = " << 0.2 << std::endl;

	ifstream fIn("field3D_64000.csv");

	std::cout << " CP = " << 0.3 << std::endl;

	string row, word;

	fIn >> row;

	std::cout << row << std::endl;

	std::cout << " CP = " << 1 << std::endl;

	_field3d.assign(nx,  vector< vector< vector< T > > >(ny, vector< vector< T > >(nz, vector< T >( 4, 0. ) ) ) );

	std::cout << " CP = " << 2 << std::endl;

	std::cout << _field3d.size() << std::endl;
	std::cout << _field3d[0].size() << std::endl;
	std::cout << _field3d[0][0].size() << std::endl;


	for(int iX = 0; iX != nx; ++iX)
	{
		for(int iY = 0; iY != ny; ++iY)
		{
			for(int iZ = 0; iZ != nz; ++iZ)
			{
				fIn >> row;

				vector<string> rawData = StringUtility::StringDelimiter(row, ",");

				StringUtility::ToTargetType<T>(rawData[1], _field3d[iX][iY][iZ][3]);
				StringUtility::ToTargetType<T>(rawData[2], _field3d[iX][iY][iZ][0]);
				StringUtility::ToTargetType<T>(rawData[3], _field3d[iX][iY][iZ][1]);
				StringUtility::ToTargetType<T>(rawData[4], _field3d[iX][iY][iZ][2]);
			}
		}
	}

	std::cout << " CP = " << 3 << std::endl;

	for(int iX = 0; iX != nx; ++iX)
	{
		for(int iY = 0; iY != ny; ++iY)
		{
			for(int iZ = 0; iZ != nz; ++iZ)
			{
				T vel[] = {_field3d[iX][iY][iZ][0], _field3d[iX][iY][iZ][1], _field3d[iX][iY][iZ][2]};
				lattice.get(iX,iY,iZ).defineRhoU(_field3d[iX][iY][iZ][3], vel);
			}
		}
	}

	std::cout << " CP = " << 4 << std::endl;

/////////////////////// Restore Field ///////////////////////


 	if(modifyBC)
 	{
		for (int iX=0; iX != nx; ++iX)
		{
			for (int iY=0; iY != ny; ++iY)
			{
				T vel[] = { 0., 0., 0.};

				lattice.get(iX,iY,   0).defineRhoU(bInput.getBoundary("Z-Left-hand-side").second, vel);
				lattice.get(iX,iY,   0).iniEquilibrium(bInput.getBoundary("Z-Left-hand-side").second, vel);

				lattice.get(iX,iY,nz-1).defineRhoU(bInput.getBoundary("Z-Right-hand-side").second, vel);
				lattice.get(iX,iY,nz-1).iniEquilibrium(bInput.getBoundary("Z-Right-hand-side").second, vel);
			}
		}
 	}

  lattice.initialize();
}

template<typename T, template<typename U> class Lattice>
void initGeometry
(
	BlockStructure3D<T, Lattice>& lattice,
	sMathLib::VectorX<int, 3> domain,
	T omega,
	Dynamics<T, Lattice>& bulkDynamics,
	OnLatticeBoundaryCondition3D<T, Lattice>& bc,
	BoundaryInput<T>& bInput
)
{
	std::cout << "LBWrapper.hh ver. 2011.5.6.1" << std::endl;
		
	const int nx = lattice.getNx();
	const int ny = lattice.getNy();
	const int nz = lattice.getNz();
	
	const int centerX = nx / 2;
	const int centerY = ny / 2;

	bool modifyBC = false;

//	T drift = 0.5;
//	T cRadius = 9.535/0.5;
//
//	std::ofstream ofstr("cRadius.txt");
//	ofstr << cRadius;
//	ofstr.close();
//
//	for(int iX = 0; iX != nx-1; ++iX)
//		for(int iY = 0; iY != ny-1; ++iY)
//			for(int iZ = 0; iZ != nz-1; ++iZ)
//			{
//				sMathLib::Vector3 rVector(iX-centerX, iY-centerY, 0.);
//
//				T r(rVector.Abs());
//
//				T dist = r - cRadius;
//
//				T phi;
//
//				if( dist < -drift )
//					phi = 0;
//				else if( dist > -drift && dist < 0)
//					phi = 1.-(dist + drift)/(drift + drift);
//				else
//					phi = 1.;
//
//				if(phi > 0)
//					lattice.get(iX, iY, iZ).defineDynamics(&instances::getNoDynamics<T, Lattice>());
//				else
//					lattice.get(iX, iY, iZ).defineDynamics(&bulkDynamics);
//			}

	lattice.defineDynamics(0,nx-1, 0,ny-1, 0,nz-1, &bulkDynamics);
 
 	if(bInput.getBoundary("Z-Left-hand-side").first  == "Pressure-Boundary")
 	{
 		bc.addPressureBoundary2N(   1,nx-2,   1,ny-2,   0,   0, omega); //set pressure B.C. on bottom
 		std::cout << "Accept the Pressure Boundary at 2N. " << std::endl;
 		modifyBC = true;
 	}
 	else if(bInput.getBoundary("Z-Left-hand-side").first  == "Velocity-Boundary")
 	{
 		bc.addVelocityBoundary2N(   1,nx-2,   1,ny-2,   0,   0, omega);
 		std::cout << "Accept the Velocity Boundary at 2N. " << std::endl;
 		modifyBC = true;
 	}

 	if(bInput.getBoundary("Z-Right-hand-side").first  == "Pressure-Boundary")
 	{
 		bc.addPressureBoundary2P(   1,nx-2,   1,ny-2,nz-1,nz-1, omega); //set pressure B.C. on top
 		std::cout << "Accept the Pressure Boundary at 2P. " << std::endl;
 		modifyBC = true;
 	}
 	else if(bInput.getBoundary("Z-Right-hand-side").first  == "Velocity-Boundary")
 	{
 		bc.addVelocityBoundary2P(   1,nx-2,   1,ny-2,nz-1,nz-1, omega);
 		std::cout << "Accept the Velocity Boundary at 2P. " << std::endl;
 		modifyBC = true;
 	}

 	//bc.addVelocityBoundary2N(   1,nx-2,   1,ny-2,   0,   0, omega);
 	//bc.addVelocityBoundary2P(   1,nx-2,   1,ny-2,nz-1,nz-1, omega);

 	bc.addVelocityBoundary0N(   0,   0,   1,ny-2,   1,nz-2, omega);
    bc.addVelocityBoundary0P(nx-1,nx-1,   1,ny-2,   1,nz-2, omega);
    //2011.10.6 modified
    bc.addVelocityBoundary1N(   1,nx-2,   0,   0,   1,nz-2, omega);
    bc.addVelocityBoundary1P(   1,nx-2,ny-1,ny-1,   1,nz-2, omega);
    
    bc.addExternalVelocityEdge0NN(   1,nx-2,   0,   0,   0,   0, omega);
    bc.addExternalVelocityEdge0NP(   1,nx-2,   0,   0,nz-1,nz-1, omega);
    bc.addExternalVelocityEdge0PN(   1,nx-2,ny-1,ny-1,   0,   0, omega);
    bc.addExternalVelocityEdge0PP(   1,nx-2,ny-1,ny-1,nz-1,nz-1, omega);


    bc.addExternalVelocityEdge1NN(   0,   0,   1,ny-2,   0,   0, omega);
    bc.addExternalVelocityEdge1NP(nx-1,nx-1,   1,ny-2,   0,   0, omega);
    bc.addExternalVelocityEdge1PN(   0,   0,   1,ny-2,nz-1,nz-1, omega);
    bc.addExternalVelocityEdge1PP(nx-1,nx-1,   1,ny-2,nz-1,nz-1, omega);

    bc.addExternalVelocityEdge2NN(   0,   0,   0,   0,   1,nz-2, omega);
    bc.addExternalVelocityEdge2NP(   0,   0,ny-1,ny-1,   1,nz-2, omega);
    bc.addExternalVelocityEdge2PN(nx-1,nx-1,   0,   0,   1,nz-2, omega);
    bc.addExternalVelocityEdge2PP(nx-1,nx-1,ny-1,ny-1,   1,nz-2, omega);

    bc.addExternalVelocityCornerNNN(   0,   0,   0, omega);
    bc.addExternalVelocityCornerNNP(   0,   0,nz-1, omega);
    bc.addExternalVelocityCornerNPN(   0,ny-1,   0, omega);
    bc.addExternalVelocityCornerNPP(   0,ny-1,nz-1, omega);
    bc.addExternalVelocityCornerPNN(nx-1,   0,   0, omega);
    bc.addExternalVelocityCornerPNP(nx-1,   0,nz-1, omega);
    bc.addExternalVelocityCornerPPN(nx-1,ny-1,   0, omega);
    bc.addExternalVelocityCornerPPP(nx-1,ny-1,nz-1, omega);

 	for (int iX=0; iX != nx; ++iX)
 	{
 		for (int iY=0; iY != ny; ++iY)
		{
 			for (int iZ=0; iZ != nz; ++iZ)
			{
				T vel[] = { 0., 0., 0. };
				lattice.get(iX,iY,iZ).defineRhoU(1., vel);
				lattice.get(iX,iY,iZ).iniEquilibrium(1., vel);
			}
		}
 	}
/////////////////////// Initialize Field ///////////////////////

 	if(modifyBC)
 	{
		for (int iX=0; iX != nx; ++iX)
		{
			for (int iY=0; iY != ny; ++iY)
			{
				T vel[] = { 0., 0., 0.};

				lattice.get(iX,iY,   0).defineRhoU(bInput.getBoundary("Z-Left-hand-side").second, vel);
				lattice.get(iX,iY,   0).iniEquilibrium(bInput.getBoundary("Z-Left-hand-side").second, vel);

				lattice.get(iX,iY,nz-1).defineRhoU(bInput.getBoundary("Z-Right-hand-side").second, vel);
				lattice.get(iX,iY,nz-1).iniEquilibrium(bInput.getBoundary("Z-Right-hand-side").second, vel);
			}
		}
 	}

 	lattice.initialize();
}

template<typename T, template<typename U> class Lattice>
void applyBoundaryPressure(BlockStructure3D<T, Lattice>& lattice, T bottomPressure)
{
	int nx = lattice.getNx();
	int ny = lattice.getNy();
	for (int iX=0; iX != nx; ++iX)
	{
		for (int iY=0; iY != ny; ++iY)
		{
			lattice.get(iX,iY,   0).defineRho(bottomPressure);
		}
	}
}


template<typename T, template<typename U> class Lattice>
void WriteTecplot(BlockStructure3D<T, Lattice> &lattice, int iT)
{
	std::cout << "Produce field data " << std::endl;	
	
	std::ostringstream ITstr; 
	
	ITstr << "field_" << iT << ".dat";
	
	std::istringstream inStream(ITstr.str());
	
	std::string filename; 
		
	inStream >> filename;
	
	DataAnalysisBase3D<T, Lattice> const &analysis = lattice.getDataAnalysis();
	
	int Nx = lattice.getNx();
	
	int Ny = lattice.getNy();
	
	int Nz = lattice.getNz();
			
	std::ofstream fluidField(filename.c_str());

	fluidField << "variables = X, Z, u, v, w, density, pressure" << endl;
	
	fluidField << "zone i = " << Nx << ", j = " << Nz << " , f = point" << endl;
		
	for(int iZ = 0; iZ != Nz; ++iZ)
	{			
		for(int iX = 0; iX != Nx; ++iX)
		{				
			fluidField << iX << " " << iZ << " ";
			
			T velocity[3] = {0., 0., 0.};
			
			velocity[0] = analysis.getVelocity().get(iX, Ny / 2, iZ)[0];
			
			velocity[1] = analysis.getVelocity().get(iX, Ny / 2, iZ)[1];
			
			velocity[2] = analysis.getVelocity().get(iX, Ny / 2, iZ)[2];
			
			fluidField << setprecision(15) << velocity[0] << " " << velocity[1] << " " << velocity[2] << " ";
			
			T pre = 0.;
				
			pre = analysis.getPressure().get(iX, Ny / 2, iZ);
				
			fluidField << setprecision(15) << pre << " " << pre << endl;
		}
	}
	
	fluidField.close();
}

template<typename T, template<typename U> class Lattice>
void WriteTecplot3DTest(BlockStructure3D<T, Lattice> &lattice, int iT)
{
	std::cout << "Produce field data " << std::endl;	
	
	std::ostringstream ITstr; 
	
	ITstr << "field_" << iT << ".dat";
	
	std::istringstream inStream(ITstr.str());
	
	std::string filename; 
		
	inStream >> filename;
	
	DataAnalysisBase3D<T, Lattice> const &analysis = lattice.getDataAnalysis();
	
	int Nx = lattice.getNx();
	
	int Ny = lattice.getNy();
	
	int Nz = lattice.getNz();
			
	std::ofstream fluidField(filename.c_str());

	fluidField << "variables = X,  Z, u, v, w, density, pressure" << endl;
	
	fluidField << "zone i = " << Nx << ", j = " << Nz << " , f = point" << endl;
		
	T fac = 1.e-3 / 1.e-4;
	
	for(int iZ = 0; iZ != Nz; ++iZ)
	{			
		for(int iY = 0; iY != Ny; ++iY)
		{			
			for(int iX = 0; iX != Nx; ++iX)
			{				
				fluidField << iX << " " << iY << " " << iZ << " ";
			
				T velocity[3] = {0., 0., 0.};
			
				velocity[0] = analysis.getVelocity().get(iX, iY, iZ)[0];
			
				velocity[1] = analysis.getVelocity().get(iX, iY, iZ)[1];
			
				velocity[2] = analysis.getVelocity().get(iX, iY, iZ)[2];
			
				fluidField << setprecision(15) << velocity[0] * fac << " " << velocity[1] * fac << " " << velocity[2] * fac << " ";
			
				T pre = 0.;
				
				pre = analysis.getPressure().get(iX, Ny / 2, iZ);
				
				fluidField << setprecision(15) << pre << " " << pre << endl;
		}
		}
	}
	
	fluidField.close();
}


template<typename T, template<typename U> class Lattice>
void WriteBinary(BlockStructure3D<T, Lattice> &lattice, int iT)
{
	std::cout << "Produce binary field data " << std::endl;		
	
	std::ostringstream ITstr; 
	ITstr << "field_" << iT << ".dat";
	std::istringstream inStream(ITstr.str());
	std::string filename; 
	inStream >> filename;
		
	DataAnalysisBase3D<T, Lattice> const &analysis = lattice.getDataAnalysis();
	
	int Nx = lattice.getNx();
	
	int Ny = lattice.getNy();
	
	int Nz = lattice.getNz();

	
	std::ofstream fOut;
	fOut.open(filename.c_str(), std::ios::out | std::ios::binary);
		
//	std::ofstream fOut_vf;
//	fOut_vf.open(filename.c_str(), std::ios::out | std::ios::binary);
		
	if(fOut.bad())
	{
		std::cerr << "Write field data error!" << std::endl;
		std::exit(0);
	}
	
	fOut.write((char*) &Nx, sizeof(int));
	fOut.write((char*) &Ny, sizeof(int));
	fOut.write((char*) &Nz, sizeof(int));

	double velocity[3];
	for(int iX = 0; iX != Nx; ++iX)
	{	
		for(int iZ = 0; iZ != Nz; ++iZ)
		{
			double velocity[3] = 
			{
				analysis.getVelocity().get(iX, Ny / 2, iZ)[0], 
				analysis.getVelocity().get(iX, Ny / 2, iZ)[1], 
				analysis.getVelocity().get(iX, Ny / 2, iZ)[2]
			};			

			fOut.write((char*)velocity, sizeof(double[3]));
			
			T pre = analysis.getPressure().get(iX, Ny / 2, iZ);
			
			fOut.write( (char*)&pre, sizeof(T) );
		}
	}

	fOut.flush();
	fOut.close();


/*
	std::ofstream fOutV, fOutP;
	std::string filenameV = "Velocity_" + filename; 
	std::string filenameP = "Pressure_" + filename; 
	fOutV.open(filenameV.c_str(), std::ios::out | std::ios::binary);
	fOutP.open(filenameP.c_str(), std::ios::out | std::ios::binary);
		
	if(fOutV.bad())
	{
		std::cerr << "Write field data (velocity) error!" << std::endl;
		std::exit(0);
	}

	if(fOutP.bad())
	{
		std::cerr << "Write field data (pressure) error!" << std::endl;
		std::exit(0);
	}
	
	fOutV.write((char*) &Nx, sizeof(int));
	fOutV.write((char*) &Ny, sizeof(int));
	fOutV.write((char*) &Nz, sizeof(int));

	fOutP.write((char*) &Nx, sizeof(int));
	fOutP.write((char*) &Ny, sizeof(int));
	fOutP.write((char*) &Nz, sizeof(int));

	const double* dTemp;
        std::vector<double> vVelocity;
        std::vector<T>      vPressure;
	int iHelfNy = Ny / 2;

	for(int iX = 0; iX != Nx; ++iX)
	{			
		for(int iZ = 0; iZ != Nz; ++iZ)
		{			
			dTemp = analysis.getVelocity().get(iX, iHelfNy, iZ);
			vVelocity.push_back(*dTemp);
			vVelocity.push_back(*(dTemp+1));
			vVelocity.push_back(*(dTemp+2));
			vPressure.push_back(analysis.getPressure().get(iX, iHelfNy, iZ));
			
		}
	}

	fOutV.write((char*)&vVelocity, sizeof(double)*3*Nx*Nz);
	fOutP.write((char*)&vPressure, sizeof(T)*Nx*Nz);


	fOutV.flush();
	fOutP.close();
	fOutV.flush();
	fOutP.close();
*/
}

template<typename T, template<typename U> class Lattice>
//void WriteVTK(BlockStructure3D<T, Lattice> &lattice, int iT, T deltaT, T deltaX, T tau)
void WriteVTK(BlockStructure3D<T, Lattice> &lattice, int iT, ibm::LBMUnit<int, T>* cu)
{
	std::cout << "Produce binary field data 3D " << std::endl;		
	
	T deltaT = cu->getDeltaT();
	T deltaX = cu->getDeltaX();
	T tau = cu->getTau();
	
	std::ostringstream osstream;
	
	osstream << "field3D_" << iT << ".vti";
		
	std::string fullName = osstream.str();
				
	DataAnalysisBase3D<T, Lattice> const& analysis = lattice.getDataAnalysis();

	T omega = 1. / tau;
	
	VtkDataWriter3D vtiOut(fullName);
	
	int nx = analysis.getVorticityNorm().getNx();
	int ny = analysis.getVorticityNorm().getNy();
	int nz = analysis.getVorticityNorm().getNz();
	vtiOut.writeHeader(0,nx-1,0,ny-1,0,nz-1,0,0,0,deltaX);
	vtiOut.startPiece(0,nx-1,0,ny-1,0,nz-1);
	//vtiOut.writeDataField(analysis.getVorticityNorm().getSerializer(IndexOrdering::backward), "vorticity", 1./deltaT, 1);
	vtiOut.writeDataField(analysis.getPressure().getSerializer(IndexOrdering::backward), "pressure", 1., 1);
	vtiOut.writeDataField(analysis.getVelocity().getSerializer(IndexOrdering::backward), "velocity", 1., 3);
	vtiOut.writeDataField(analysis.getVorticity().getSerializer(IndexOrdering::backward), "vorticity", 1., 3);
  
	T scale = -omega / (T)2 * Lattice<T>::invCs2;//1/scale
	vtiOut.writeDataField(analysis.getStrainRateFromStress().getSerializer(IndexOrdering::backward), "stress", 1., 6);
	vtiOut.endPiece();
	vtiOut.writeFooter();
  
	analysis.reset();
  
	
	//DataSerializer<T> const& serializer = analysis.getVelocity().getSerializer(IndexOrdering::forward);
	//saveAsciiData<T>(analysis.getPressure(), filename);
	/*
	std::ofstream fOut(filename.c_str());
	serializer2ostr<T>(serializer, &fOut, false);
	*/
/*
	int Nx = lattice.getNx();
	
	int Ny = lattice.getNy();
	
	int Nz = lattice.getNz();

	
	std::ofstream fOut;
	fOut.open(filename.c_str(), std::ios::out | std::ios::binary);
		
	std::ofstream fOut_vf;
	fOut_vf.open(filename.c_str(), std::ios::out | std::ios::binary);
		
	if(fOut.bad())
	{
		std::cerr << "Write field data error!" << std::endl;
		std::exit(0);
	}
	
	// 1. write size of field and particle number
	fOut.write((char*) &Nx, sizeof(int));
	fOut.write((char*) &Ny, sizeof(int));
	fOut.write((char*) &Nz, sizeof(int));

	double velocity[3];
	for(int iX = 0; iX != Nx; ++iX)
	{
		for(int iY = 0; iY != Ny; ++iY)
		{
			for(int iZ = 0; iZ != Nz; ++iZ)
			{
				double velocity[3] 
				
				
				// 2. write ux, uy, uz
				fOut.write((char*)velocity, sizeof(double[3]));
			
				T pre = analysis.getPressure().get(iX, iY, iZ);
			
				// 3. write p
				fOut.write( (char*)&pre, sizeof(T) );				
			}
		}
	}

	fOut.flush();
	fOut.close();
*/

/*
	std::ofstream fOutV, fOutP;
	std::string filenameV = "Velocity_" + filename; 
	std::string filenameP = "Pressure_" + filename; 
	fOutV.open(filenameV.c_str(), std::ios::out | std::ios::binary);
	fOutP.open(filenameP.c_str(), std::ios::out | std::ios::binary);
		
	if(fOutV.bad())
	{
		std::cerr << "Write field data (velocity) error!" << std::endl;
		std::exit(0);
	}

	if(fOutP.bad())
	{
		std::cerr << "Write field data (pressure) error!" << std::endl;
		std::exit(0);
	}
	
	fOutV.write((char*) &Nx, sizeof(int));
	fOutV.write((char*) &Ny, sizeof(int));
	fOutV.write((char*) &Nz, sizeof(int));

	fOutP.write((char*) &Nx, sizeof(int));
	fOutP.write((char*) &Ny, sizeof(int));
	fOutP.write((char*) &Nz, sizeof(int));

	const double* dTemp;
        std::vector<double> vVelocity;
        std::vector<T>      vPressure;
	int iHelfNy = Ny / 2;

	for(int iX = 0; iX != Nx; ++iX)
	{			
		for(int iZ = 0; iZ != Nz; ++iZ)
		{			
			dTemp = analysis.getVelocity().get(iX, iHelfNy, iZ);
			vVelocity.push_back(*dTemp);
			vVelocity.push_back(*(dTemp+1));
			vVelocity.push_back(*(dTemp+2));
			vPressure.push_back(analysis.getPressure().get(iX, iHelfNy, iZ));
			
		}
	}

	fOutV.write((char*)&vVelocity, sizeof(double)*3*Nx*Nz);
	fOutP.write((char*)&vPressure, sizeof(T)*Nx*Nz);


	fOutV.flush();
	fOutP.close();
	fOutV.flush();
	fOutP.close();
*/
}

#endif //_LBWRAPPER_HH_
